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SUMMARY 


The  present  study  addresses  the  aerodynamics  of  a  lobed  body  immersed  in  a  Mach  6 
hypersonic  flow  field  at  sea  level.  As  a  first  study  of  this  type,  the  shape  of  the  body  is  held 
fixed,  and  the  flow  field  is  resolved  by  applying  state-of-the-art  large  eddy  simulation  techniques 
in  conjunction  with  a  hybrid  shock-turbulence  capturing  algorithm.  Air  is  treated  as  a  mixture  of 
nitrogen  and  oxygen,  and  the  governing  equations  are  closed  by  a  modern  compressible 
turbulence  closure  term.  Pressure  is  determined  by  using  the  thermally  perfect  gas  equation  of 
state  applied  to  each  species.  The  distribution  of  temperature  is  determined  on  the  body  surface 
as  well  as  temperature  gradients  based  upon  adiabatic  wall  boundary  conditions.  The  structure  of 
the  flow  field  is  examined  as  is  the  time  required  for  stationarity.  Both  turbulent  statistics  and 
spectra  are  determined  at  points  on  the  surface  of  the  body.  Times  series  analyses  are  performed 
for  aerodynamics  forces  encountered  by  the  body. 
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1  INTRODUCTION 


The  past  two  decades  have  been  witness  to  major  advances  in  computing  power,  particularly 
in  the  area  of  parallel  computing.  For  instance,  in  1995,  the  author’s  doctoral  research  was 
performed  on  a  CRAY  Y-MP  with  vector  registers,  512  MB  of  memory  and  a  clock  speed  of  167 
MHz.  The  system  had  no  more  than  eight  vector  processors  and  had  less  capability  than  the 
ThinkMate  8-Way  system  (32  cores)  that  resides  by  my  desk.  Now,  we  have  high  performance 
computers  (HPC)  platforms  such  as  RAPTOR,  a  CRAY  XE6  with  40,000  cores  and  2GB 
memory  per  node.  There  are  32  cores  per  node  making  for  a  staggering  improvement  in 
computational  capability.  This  windfall  in  computer  power  constitutes  a  boon  for  the  field  of 
computational  physics.  Modest  improvements  in  numerical  algorithms,  particularly  those  for 
solving  partial  differential  equations  (PDEs),  can  now  be  fully  implemented  and  exercised  with 
fewer  restrictions  imposed  by  memory  size  or  execution  speed.  Fields  such  as  Computational 
Electromagnetics  (CEM),  Computational  Structural  Mechanics  (CSM)  and  Computational  Fluid 
Dynamics  (CFD)  can  benefit  from  the  use  of  research-grade  numerical  schemes  and  state-of-the- 
art  material  models.  In  some  cases,  computational  tasks  requiring  several  weeks  of  computer 
time  (just  five  years  ago)  can  now  be  completed  in  a  few  hours.  We  do  not  exaggerate  when 
today  we  claim  that  many  chemical  dynamics  and  material  modeling  problems  can  be  solved  via 
Quantum  Chemistry  algorithms  on  HPC  systems.  Budgetary  considerations  aside,  it  is  a 
fortuitous  time  to  work  in  the  field  of  Computational  Physics. 

From  an  egocentric  standpoint,  the  field  of  CFD  is  undergoing  a  transition.  Aside  from  its 
earlier  focus  on  solving  problems  in  aerodynamics,  CFD  is  now  working  to  incorporate 
additional  physics  within  its  algorithmic  bag  of  tricks.  Some  of  the  more  recent  and  important 
work  enables  both  shock  waves  and  turbulent  eddies  to  be  captured  without  mutual 
degradation.  [1,2]  State-of-the-art  dynamic  methods  for  Farge  Eddy  Simulation  (FES)  of 
turbulence  are  also  carried  into  the  realm  of  shock  physics  as  a  unified  algorithm.  This  advance 
presents  a  powerful  toolset  for  capturing  the  physics  of  extreme  environments  such  in  the  violent 
interior  of  stars  or  in  the  core  of  manmade  explosions. [3,4]  The  current  research  field,  regarded 
hereafter  as  Multiphase  Physics  (MP),  concentrates  on  simulating  matter  fields  of  just  level  of 
complexity  where  multiple  phases  of  matter  may  be  involved  along  with  chemistry  and  extreme 
pressures  and  temperatures.  What  is  extreme?  From  the  standpoint  of  some  of  our  recent  work, 
extreme  may  be  defined  as  temperatures  above  500CFK  or  pressures  cast  in  the  Giga-Pascal 
range.  Naturally,  MP  problems  are  quite  realistic  in  that  they  may  contain  solid  particles  or  liquid 
droplets  that  combust  or  react  chemically  with  a  surrounding  gaseous,  turbulent  flow  field. [5] 
From  our  point  of  view,  a  properly  formulated  set  of  numerical  algorithms  contains  all  of  the 
requisite  physics  for  the  problem;  these  algorithms,  in  turn,  autonomously  resolve  the  problem 
physics  based  upon  only  the  input  parameters  such  as  geometry  and  initial  conditions  (including 
the  relevant  chemical  species).  It  is  out  desire  to  minimize  or  even  eliminate  any  processes  such 
as  “tuning”  of  “calibration”  of  the  computer  algorithms  undertaken  on  the  part  of  the  user.  The 
problem  under  study  below  exemplifies  this  idea. 

As  we  alluded  to  earlier,  the  physics  of  extreme  environments  is  very  interesting.  The 
hypersonic  flow  field  created  around  a  solid  body  under  sea  level  flight  conditions  is  an  example 
of  just  such  an  environment.  Although  hypersonic  flow  was  of  great  interest  in  the  final  quarter 
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of  the  twentieth  century,  the  focus  of  the  research  was  largely  confined  to  high  altitude  flight. 
This  effort  was  motivated  by  the  space  program.  We  are  concerned  with  the  behavior  of 
hypersonic  flow  in  the  thick  atmosphere  where  gases  such  as  oxygen  and  nitrogen  have  higher 
partial  pressures.  In  combination  with  trace  gases,  these  species  are  likely  to  react  chemically 
near  the  surface  of  an  immersed  body.  As  one  may  expect,  temperatures  will  soar  in  the  body’s 
boundary  layer  as  viscous  forces  slow  the  flow  stream  to  satisfy  the  no-slip  boundary  condition. 
The  variation  in  temperature  can  vary  over  a  wide  range,  so  our  equation  of  state  (EOS)  must 
correctly  capture  changes  in  the  gaseous  specific  heats.  It  is  worth  noting  that  the  older 
calorically  perfect  gas  model  over  predicts  temperatures.  During  the  period  of  time 
encompassing  the  1960s  through  the  1980s,  hypersonic  flow  and  high  temperature  gas  dynamics 
were  intensively  studied  because  of  the  design  needs  of  reentry  vehicles  and  concept  designs 
such  as  the  National  Aerospace  Plane.  [6]  In  deference  to  supersonic  flow,  there  is  no  strict  point 
of  demarcation  in  terms  of  speed  heralding  the  hypersonic  regime.  For  aerospace  vehicles 
commonly  studied,  the  onset  of  hypersonic  flow  physics  is  encountered  in  the  range  of  Mach  3  to 
Mach  5,  or  even  higher.  The  physics  associated  with  this  regime  is  very  different  that  that 
associated  with  slower  speed  flow.  Hypersonic  flow  is  characterized  by  the  presence  of  thin 
shock  layers  that  lie  between  the  body  and  the  shock  wave. [6]  At  higher  Mach  numbers,  this 
layer  becomes  even  thinner,  and  it  may  even  merge  with  the  boundary  layer.  The  presence  of  an 
entropy  layer  is  also  characteristic  of  many  hypersonic  bluff  body  flow  fields.  Bluff  (or  blunted) 
bodies  cause  the  formation  of  curved  shock  waves  that  stand  some  distance  away  from  the  body. 
The  strength  of  the  shock  wave  and  the  attendant  entropy  jump  changes  along  with  shock 
curvature.  It  is  strongest  at  the  “nose”  and  weakens  as  the  shock  wave  straightens  downstream. 
This  effect  causes  a  flow  with  strong  entropy  and  vorticity  gradients  to  wash  down  the  body 
thickening  the  boundary  layer.  [6]  High  temperature  gas  dynamics  are  often  encountered  in 
hypersonic  flow  fields.  In  order  to  satisfy  the  no-slip  boundary  condition,  the  flow  must  reduce 
from  hypersonic  speed  outside  of  the  boundary  layer  to  zero  at  the  body  surface.  This  situation 
stands  in  stark  contrast  to  earlier  models  of  hypersonic  flow  that  allow  slip  (or  tangential  flow)  at 
the  body  surface.  The  mechanism  decreasing  the  flow  speed  is  viscous  dissipation  or  friction. 
The  kinetic  energy  of  the  flow  is  transformed  into  the  internal  energy  of  gas  molecules  in  the 
boundary  layer.  As  a  result,  boundary  layer  temperature  increases  dramatically.  [7]  In  slower 
speed  flow  fields,  this  energy  is  absorbed  mostly  in  molecular  translational  and  rotational  modes, 
but  for  hypersonic  speeds,  energy  may  be  absorbed  in  vibrational  modes. [8]  The  existence  of  the 
latter  mode  tends  to  moderate  the  rise  in  temperature  and  motivates  the  use  of  more  advanced 
EOS  models  for  gas  mixtures. 

Over  more  complex  body  shapes,  the  attendant  high  Reynolds  number  leads  us  to  expect  the 
formation  of  pockets  of  turbulent  flow.  It  is  important  that  our  physics  simulation  capture  this 
behavior  since  turbulence  has  a  potent  effect  over  chemical  reactions.  In  turn,  chemical  reactions 
affect  temperature  at  the  body  surface.  All  of  these  effects  combine  to  shape  the  aerodynamic 
forces  (or  loads)  exerted  on  the  body.  The  loads,  of  course,  exert  a  great  deal  of  influence  over 
the  body’s  flight  stability.  A  goal  of  this  report  is  to  investigate  the  evolution  of  turbulent  flow  at 
hypersonic  speeds  under  sea  level  flight  conditions.  It  is  desirable  to  consider  an  immersed  body 
with  enough  complexity  to  support  turbulent  flow  while  retaining  the  ease  of  single  block  grid 
generation.  As  is  shown  in  Figure  1,  our  test  body  has  a  complex  shape  characterized  by  a  blunt 
nose  and  an  elongated  body  distending  into  fin-shaped  lobes.  The  body  is  immersed  in  a  Mach  6 
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flow  field  with  angle  of  attack  and  side  slip  angle  both  fixed  at  5°.  The  grid  for  this  body  is 
generated  through  the  use  of  algebraic  techniques;  the  specifics  are  discussed  later  in  this  report. 


Figure  1.  Single-Block  Test  Geometry 


State-of-the-art  LES  techniques  are  used  to  simulate  flow  around  the  body  at  the 
aforementioned  flight  conditions.  These  methods  allow  the  viscous,  turbulent  flow  field  to  be 
resolved  with  high  accuracy.  Adiabatic  wall  conditions  are  enforced  at  the  solid  surface  for  the 
resolution  of  temperature.  In  the  course  of  the  simulation,  the  time  required  to  achieve  a 
stationary  flow  field  is  determined.  Flow  properties  vary  in  a  cyclical  manner  indicative  of  a 
stationary  flow.  Pressure,  sub  grid  kinetic  energy,  temperature  and  temperature  gradient 
distributions  are  also  determined  on  the  body  surface.  Velocity  correlations,  turbulent  kinetic 
energy  spectra  and  probability  distribution  functions  for  fluctuating  velocity  components  are 
determined  at  points  near  the  body  surface. 

The  remainder  of  this  report  is  organized  as  follows.  Section  2  contains  an  exposition  of  the 
theory  behind  our  LES  methodology  including  current  research  in  subgrid  models.  Section  3 
contains  brief  descriptions  of  the  numerical  techniques  used  in  LESLIE3D  as  well  as  a 
discussion  of  our  grid  generation  technique.  The  statistical  analysis  techniques  employed  are  also 
briefly  presented.  Some  of  the  specifics  associated  with  the  set-up  of  our  simulation  are 
presented  in  Section  4.  Section  5  contains  the  results  of  our  simulations  at  Mach  6  along  with  the 
analyses  produced  by  statistical  post-processing.  Section  6  of  the  report  is  dedicated  toward 
conclusions  along  with  a  wrap-up  of  our  findings. 

2  THEORY 

The  numerical  flow  field  solutions  described  below  have  been  generated  by  using  the  Large 
Eddy  Simulation  with  Linear  Eddy  modeling  in  3  Dimensions  (LESLIE3D)  computer  program. 
LESLIE3D  is  a  state-of-the-art  MP  research  tool  developed  by  Suresh  Menon  at  the  Georgia 
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Institute  of  Technology.  [9]  This  computer  program  has  a  core  capability  for  dynamic  large  eddy 
simulation  (LES).  That  is  to  say,  flow  field  features  existing  at  the  scale  of  the  grid  or  larger  are 
simulated  by  numerical  solution  of  the  conservation  equations.  Flow  features  existing  at  the 
subgrid  scales  (scales  smaller  than  the  mesh  size)  are  modeled. [10]  The  subgrid  effects  are 
represented  in  the  conservation  equations  via  closure  terms.  The  theory  behind  LESLIE3D  is 
described  in  the  following  sections. 


2.1  Filtered  Governing  Equations 


The  gas  phase  conservation  laws  used  in  LESLIE3D  consist  of  the  compressible  Navier- 
Stokes  equations  cast  in  three  dimensions.  For  real  turbulent  flow  fields  (in  spite  of  today’s 
computer  resources),  we  cannot  solve  these  equations  via  direct  numerical  simulation  (DNS)  for 
Reynolds  numbers  exceeding  8000.  [11]  One  may  recall  that  hypersonic  flow  fields  are 
characterized  by  Reynolds  number  of  one  million  or  more.  Our  LES  approach  first  involves 
spatially  filtering  the  governing  equations  at  the  grid  scale  in  order  to  produce  a  set  of  turbulence 
closure  terms  that  are  later  modeled.  [12]  The  Navier-Stokes  equations  are  filtered  in  order  to 
permit  the  large  spatially  dependent  scales  for  fluid  motion  to  be  separated  from  the  small 
universal  scales. [13]  In  the  equations  below,  the  overbar  indicates  a  spatially  filtered  quantity 
while  the  tilde  (~)  indicates  a  mass  averaged  property,  i.e., 


P 


(1) 


The  filtered  Navier-Stokes  equations  may  be  written 


dp  |  dpu,  _  Q 

dt  dxi 


(2) 


dpu, 

dt 


+  -^—\OU.U.  +  PS; 

dXj  '  J 


sgs  1  _ 


r,..  +  rr;"] 


0 


(3) 


dpE 

dt 


d  _  ~  —  _  _ 

+  ~ — [ (pE  +  P)  Uj  +  qt 
dx, 


Ujhi 


+  + of' ]  =  0 


(4) 


dpYk 

dt 


d 

+  —lpYk(p 

ox 


+v!,k)+Y,7  +  0-is]  =  ak 


(5) 


where  k  =  1, . . .,  Ns ,  Ns  being  the  number  of  gaseous  chemical  species  involved  in  the  calculation. 
Equations  (2)  through  (5)  are  mass,  momentum,  energy  and  species  conservation  equations, 
respectively.  Symbols  p  ,  ui  ,  P  ,  Yk  and  q;  are  the  mixture  gas  density,  velocity  component  in 

Cartesian  component  direction  xf,  absolute  pressure,  mass  fraction  for  the  k'h  species  and  heat 
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flux  component  in  direction  x;,  respectively.  Symbol  rij  contains  the  components  of  the 
Cartesian  shear  stress  tensor  while  E  represents  the  total  energy  per  unit  volume,  i.e., 


Tn=M 


f  duj  duj  ^ 
Kdxj  dxt  j 


+A:l 

Sx,  11 


(6) 


E  =  e  +  —  u,  ul 


(7) 


The  summation  convention  is  applied  to  the  index  /  in  (6)  and  (7).  The  Fickian  diffusion 
velocities  are  given  by  Vt  k .  The  subgrid  stress  tensor  tEs  ,  subgrid  total  enthalpy  flux  H’“ , 

subgrid  convective  species  flux  YtJs ,  subgrid  viscous  work  erf®5  and  subgrid  species  diffusive 
work  6-sks  are  subgrid  quantities  produced  by  the  filtering  process.  Their  exact  forms  may  be 
written  as 


Tij=P(<  Uiuj>-uiuj) 

(8) 

=  p(<E  ut  >  —E  ui  )  +  UjP  —  P  Uj 

(9) 

(10) 

<TF=-(.UjTii-UjTii) 

(11) 

er  =  p«VikYk>-VikYk) 

(12) 

The  filtered  species  diffusion  velocities  are  given  by 


v.  =-^ay‘ 


yk  ax,. 


(13) 


where  Dk  is  the  diffusion  coefficient  for  species  k.  In  equations  (8)  through  (12),  the  arrowhead 

brackets  “<  >”  denote  mass  averaging  as  does  the  notation.  These  equations  contain  averages 
of  products  between  variables.  The  repeated  index  in  (12)  does  not  imply  summation.  Not  being 
known  a  priori ,  the  subgrid  terms  must  be  treated  as  variables.  The  result  is  that  the  system  of 
filtered  governing  equations  now  has  more  variables  than  equations  and  cannot  directly  be  closed 
for  solution.  Closure  can  only  be  accomplished  by  modeling  these  subgrid  quantities.  [10] 

Thermochemical  and  thermophysical  behavior  for  the  system  must  also  be  described. 
Chemical  reactions  occurring  between  species  are  governed  by  a  reaction  mechanism  specified 
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by  the  user.  The  rate  of  mass  transfer  of  between  species  is  given  by  cok ,  the  filtered  reaction  rate 

term.  [12]  Pressure  and  temperature  in  the  gas  mixture  are  determined  through  the  use  of  a 
filtered  perfect  gas  equation  of  state,  i.e., 


P  =  p(RT+Tsgs) 


with  Rl,  equal  to  the  universal  gas  constant,  and 


k= 1 


*'  /  Rv  A 


k= 1 


(14) 


(15) 


MWk  is  the  molecular  weight  of  species  k  .  The  sum  of  subgrid  species-temperature  correlations 
is  given  by  Tsgs,  i.e., 


r,;  f <YkT>-YkT 

h  MWk 

The  filtered  mixture  internal  energy  may  be  expressed  as 

?  =  +  E  ?.  ]cv.t  (T), rfr  +  f;  El 

k= 1  k= 1  0  *=1 


(16) 


(17) 


In  equation  (9),  Cv,  t  is  the  constant  volume  specific  heat  for  species  k ,  and  Ekgs  is  the  subgrid 
internal  energy  for  species  k  ,  i.e., 


Er=<YkeAT)>-YkeAT) 


(18) 


The  heat  flux  term  may  be  written  as 


prr  Ns  Ns 

q,=-K(f)  —  +  pYJhk(f)YlVil+t.  «3' 

OX;  k=1  k=1 


(19) 


where  k  is  the  local  coefficient  of  thermal  conductivity.  [9]  The  subgrid  species  heat  flux  may  be 
expressed  (with  no  summation  over  k)  as 


qT  =  P«hkYkVik>-hkYkVik ) 


(20) 


Note  that  repeated  indices  in  (18)  and  (20)  do  not  imply  summation.  Also  the  filtered  sensible 
species  enthalpy  is 
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(21) 


f 

K=yth “i+y,Jc„(r)rfr 

0 

The  specific  heat  at  constant  pressure  for  species  k  is  written  as 

Cp,k  =  CV'k  + 


(22) 


2.2  Modeling  Correlations  to  Obtain  Closure 


The  subgrid  terms,  denoted  by  superscript  “ sgs ”  cause  difficulties  in  closing  the  system  of 
governing  equations.  LES  fixes  this  problem  by  modeling  these  terms  based  upon  the  resolved 
flow  properties.  Specifically,  we  neglect  Tsg\  9-gs  ,  q*gs  and  Eskgs  since  these  terms  tend  to  be 

small.  [12]  The  crux  of  the  modeling  process  lies  in  the  determination  of  the  subgrid  stress  tensor 
based  upon  the  sub  grid  kinetic  energy,  i.e., 


r 


t;*'  =  -2pv, 


S ij 


'A  2 

+  —ksgs5i . 
J  3  ,J 


(23) 


where  the  evolution  equation  for  subgrid  kinetic  energy  is  given  by 


d  d  d 

^(pkn+^ipu^n  =  f- 

Ot  OX:  OX: 


dksgs  pv,R8T 

(y OV.  +  p) - +  — - 

dxl  Pr;  dXj 

* 


l  +  aAM^)2 


pSks 

Dsgs 

V  Uk  J 


\  2  A 


r 


Bu¬ 


rn 


(yts^)3/23 


T-gS — -  + 

v  7  dxt  A  , 


The  subgrid  kinetic  energy  is  mathematically  defined  as 

ksgs  =  ^(<  ujui  >  -u,u, ) 


(25) 


where  the  summation  convention  is  in  effect  for  the  index  1  .[9]  Equation  (24)  is  integrated  along 
with  the  conservation  equations,  and  as  a  model  equation,  relies  on  a  set  of  parameters  such  as 
the  eddy  viscosity  v,  given  by 


vt=KvJk^A  (26) 

where  Kr  is  a  coefficient  set  by  the  Locally  Dynamic  subgrid  Kinetic  energy  Model  (LDKM).[9] 
The  dissipation  term  for  ksgs  is 
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with  Ks  as  the  LDKM  dissipation  coefficient  (analogous  to  Ks );  A  is  the  local  grid  scale 
measurement.  For  most  problems,  Krand  K  arc  set  to  0.067  and  0.931,  respectively.  Other 
parameters  introduced  in  (24)  are  ju  the  dynamic  viscosity,  Pr(  the  turbulent  Prandtl  number, 
Mstgs  the  turbulent  Mach  number  based  upon  ksgs ,  and  a  d  is  a  pressure-dilatational  scaling 
coefficient.  Note  also  that  (23)  and  (24)  are  tightly  coupled  through  the  presence  of  both  ksgs  and 


The  theory  that  underlies  LDKM  is  quite  complex,  and  the  details  of  implementing  the  above 
equations  (along  with  others  not  shown)  are  beyond  the  scope  of  this  report.  This  statement  is 
especially  true  when  we  employ  the  algorithms  developed  for  dynamically  updating  K,,  and  K  .; . 
Still,  it  is  worthwhile  to  discuss  the  motivation  behind  this  model.  In  the  first  place,  one  may  take 
note  of  the  number  of  coefficients  existing  in  (24),  (26)  and  (27).  The  proliferation  of 
coefficients  indicates  that  LDKM  is  indeed  a  model,  but  these  coefficients  are  not  tunable. 
Instead,  they  are  in  most  cases  based  upon  knowledge  of  the  dynamic  behavior  of  turbulent 
compressible  flow  fields.  That  is  to  say,  these  coefficients  are  set  either  autonomously  by  LDKM 
operating  in  dynamic  mode  or  remain  fixed  as  described  above.  [12]  Also,  LDKM  diverges  from 
most  contemporary  dynamic  LES  models  through  the  presence  of  the  term  marked  by  an  (*) 
asterisk.  This  term  represents  dilatational  effects  within  the  flow  field  due  to  pressure.  It  is 
through  this  term  that  (24)  becomes  representative  of  highly  compressible  flow  fields.  This 
equation  stands  as  an  example  of  state-of-the-art  research  in  this  field.  When  operating  LDKM  in 
dynamic  mode,  the  subgrid  scale  properties  are  made  to  behave  in  an  inertially  correct  manner 
with  respect  to  the  resolved  flow  field. 

The  remainder  of  the  subgrid  terms  are  not  as  difficult  to  model  although  the  subgrid  kinetic 
energy  is  still  relied  upon  (at  least  in  part).  H.gs  and  cr.gs  are  modeled  as  one  by  using  eddy 
viscosity  closure  form,  i.e., 

r)ksgs  OV  C  r)T 

Hsgs  +  of s  =  ~(p  v,  +  p)  — - — — — — i-  u ;  rf  (28) 

'  '  '  ax,.  Pr,  ax,.  7  " 

The  subgrid  diffusion  of  species  mass  fractions  is  modeled  by  using  an  eddy  diffusivity  form, 
i.e., 


y  sgs  _  _  P  K 

Sct  ax,. 
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where  Sc,  is  the  turbulent  Schmidt  number.  The  dynamic  viscosity  jj  is  computed  from 
Sutherland’s  law.  That  is, 


for  a  reference  viscosity  go- 


f  T  A3/2 
VTo  J 


T0+S 
T  +  S 


(30) 


3  NUMERICAL  METHODS 

In  the  preceding  section,  the  filtered  Navier-Stokes  equations  are  presented  along  with  the 
equations  for  LDKM.  These  equations  are  well  suited  for  describing  turbulent  flow  fields  in  any 
flight  regime  (subsonic,  supersonic  and  hypersonic).  In  this  part  of  the  report,  we  discuss  the 
numerical  algorithms  associated  with  solving  these  equations.  Our  grid  generation  method  is  also 
described.  The  hypersonic  flow  problem  of  interest  is  characterized  by  the  presence  of  strong 
shock  waves  and  high  temperature  gas  dynamics.  Historically,  it  has  proven  difficult  for 
computer  codes  to  simultaneously  capture  both  shock  waves  and  subtle  turbulent  flow  field 
fluctuations.  This  difficulty  is  caused  by  differences  between  high  accuracy  space  schemes  and 
schemes  designed  to  capture  strong  discontinuities.  Generally,  numerical  dissipation  is  used  to 
dampen  oscillations  around  shock  waves,  yet  this  same  dissipation  washes  away  the  weaker 
vortical  motions  associated  with  turbulence.  On  the  other  hand,  high  order  centered  difference 
stencils  are  usually  required  to  resolve  turbulence  eddies,  but  these  stencils  cause  severe  spurious 
oscillations  around  shock  waves.  It  is  also  important  to  realize  that  shock-capturing  algorithms 
must  be  able  to  contend  with  complex  equations  of  state.  For  instance,  we  prefer  to  use  the 
thermally  perfect  gas  equation  of  state  in  the  hypersonic  flight  regime  in  order  to  accurately 
resolve  temperature.  To  satisfy  these  requirements,  our  efforts  have  produced  a  two-part  scheme. 
For  shock-capturing,  we  apply  what  is  referred  to  as  the  Harten-Lax-van  Leer-Contact/Einfeldt 
(HLLC/E)  approximate  Riemann  solver. [14, 15]  In  smooth  regions  of  the  flow  field,  we  apply  a 
second/fourth  order  MacCormack  solver  used  in  the  earlier  incarnations  of  LESLIE3D.  A 
switching  algorithm  allows  the  computer  program  to  autonomously  select  the  appropriate  space 
scheme  based  upon  local  properties  throughout  the  flow  field.  Below,  we  present  brief 
descriptions  of  the  spatial  integration  schemes. 


3.1  The  HLLC/E  Approximate  Riemann  Solver 

HLLC/E  is  composed  of  HLLE  (HLL-Einfeldt),  a  scheme  that  does  not  preserve  the  contact 
wave,  and  HLLC  (HLL-Contact),  a  scheme  that  does  capture  the  contact  wave.  In  this  context, 
the  contact  wave  is,  in  fact,  a  contact  discontinuity,  a  common  flow  feature  exhibited  by  the 
shock  tube  problem  as  well  as  by  strong  blast  waves.  These  numerical  schemes  are  constructed 
on  a  one-dimensional  coordinate  oriented  normal  to  a  finite  volume  interface.  We  can  orient  this 
coordinate  across  any  of  a  finite  volume  cell’s  interfaces,  so  the  technique  is  readily  usable  in 
three  dimensions.  The  HLLC/E  solution  is  calculated  for  a  vector  of  conserved  variables 


U  =(p,u,v,w, pE , pksgs, pYu..., pYNs ) 


(31) 


the  associated  flux  vector  is 
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F  =  (pq, puq  +  Pnx,pvq  +  Pny , pwq  +  Pnz , (pE  +  P)q, pksgsq , pYAq,.. ., pYN  q)  (32) 


where  q  =  un+vnv  +wn  ,  and  f=ni  +n„  /  +  n  k  is  the  unit  surface  normal  vector  at  the 

J-  x  y  z  '  x  y  */  z 

interface. [12]  Based  upon  these  definitions,  the  HLLE  numerical  flux  vector  can  be  written  as 


t?HLLE 

j+1/2 


F\  0  <  S1 

<  F\  SL  <0<  SR 
Fr,  Sr<  0 


(33) 


where 


SRFL-SLFR  +SLSR(UR-UL) 

F  = - FT? -  (34) 

One  may  note  from  the  form  of  (33),  that  the  HLLE  method  recognizes  a  shock  discontinuity, 
but  not  a  contact  surface.  The  Einfeldt  wave  speed  estimates  are 

SL  =  min{qL  -cL,q  —  c)  (35) 


SR  =ma x(qR  +cR  ,q  +  c) 


(36) 


In  (35)  and  (36),  c  is  the  speed  of  sound,  and  the  notation  (  “  )  indicates  Roe  averaging. [16]  This 
numerical  flux  works  quite  well  for  hypersonic  flow  and  is  not  susceptible  to  contact  instabilities 
that  can  occur  when  strong  bow  shocks  are  aligned  with  the  grid.  Unfortunately,  because  of 
dissipative  effects  near  the  contact  wave,  it  cannot  be  used  for  extended  regions  of  the  flow  field 
where  turbulence  may  be  present.  [12]  It  is  also  best  not  to  apply  this  scheme  for  moving  shock 
waves  or  obliquely  oriented  shocks.  Under  these  circumstances,  the  HLLC  numerical  flux  shows 
improved  performance.  The  numerical  flux  for  the  HLLC  scheme  has  a  different  form,  i.e., 


rHLLC 

U+1/2 


FL 

fl*  =  Fl  +  Sl(Ul* -Ul), 
fr*  =  Fr  +  Sr(Ur*  -Ur), 
FR 


0  <SL 
SL  <0  <  s* 
SL*  <0<  SR 
SR<  0 


(37) 


In  (37),  intermediate  states  UL>  and  UR'  are  introduced  in  order  to  model  the  presence  of  the 
contact  discontinuity.  [14]  These  states  may  be  written  as 

UL  =  aLUL  +  (0,pLcoLnx,pLcoLny,pLCQLnz,i/fL  ,0,...)  (38) 
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with 


C*  L 

PL=  L  q* 

SL-S 

(39) 

+ 

II 

a 

(40) 

coL=-/3L(qL-SL ) 

(41) 

L  PS  —  PLqL 

V  = - 7 - P— 

SL-S 

(42) 

and  the  intermediate  wave  speed 

j~)R  j~)L  .  _L^L/rrL  \  ^ rr  /?  /?  \ 

o*_-P  -p  +/0  <?  (5  -g  )-/?  q  (5  -g  ) 
p  (S  -q  )-p  (S  -q  ) 


The  state  [7s  is  formed  by  replacing  R  for  L  in  (38)  through  (43).  In  LESLIE3D,  both  of  these 
numerical  flux  formulations  are  used  to  form  the  HLLC/E  scheme.  We  state  the  result  only,  i.e., 


77  HLLC  /  E 

U+1/2 


I  Fitv2  ’  (dpj  <  0  and  du  j  <  0)  or  (dp  k  <  0  and  du  k  <  0) 
1  Fitv2  >  otherwise 


(44) 


This  method  requires  shock  detection  in  directions  transverse  to  the  direction  of  the  interface 
normal,  so 


I  Pi,j+-\,k  ^j-u  I  _  1 

u)  3 


(45) 


du.j  =  uij+ u 


-i,t 


(46) 


and  so  forth. 


A  close  inspection  of  the  numerical  flux  formulas  above  reveals  that  the  computation  of  the 
numerical  flux  first  requires  the  left  (L)  and  right  ( R )  states  be  identified  at  (i  + 1/2),  an 
interface  location.  These  “upwind”  states  are  reconstructed  from  the  data  stored  in  the  finite 
volume  cells.  This  process  is  accomplished  with  the  use  of  Monotone  Upwind  Scheme  for 
Conservation  Laws  (MUSCL)  interpolation. [17]  High  order  interpolation  also  requires  the  use  of 
a  nonlinear  limiter  with  “flattening”  in  order  to  maintain  data  monotonicity.  [18]  The  details  of 
these  procedures  are  omitted  from  this  work. 
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3.2  The  MacCormack  Solver 

For  the  “hybrid”  shock-turbulence  capturing  solver,  it  is  necessary  to  switch  between  the 
upwind  scheme  and  a  centralized  scheme  based  upon  MacCormack ’s  method  when  calculating 
the  numerical  flux  at  interfaces.  The  finite  difference  version  of  this  scheme  uses  forward  and 
backward  differences  alternatively  to  remove  bias  from  numerical  error.  For  the  finite  volume 
scheme,  we  choose  forward  and  backward  upwind  variables  alternatively  for  the  flux 
computation,  i.e., 


= m iv2) 


(47) 


where 


UUi2=UM  ; 


U;+v2=Ui 


(48) 


It  can  be  shown  that  this  scheme  is  stable  and  second  order  accurate  in  space.  LESLIE3D  also 
offers  a  MacCormack  space  scheme  accurate  to  the  fourth  order.  The  hybrid  method  is  required 
to  switch  back  and  forth  between  the  upwind  and  MacCormack  solvers.  The  “switch”  is  based 
upon  the  “smoothness”  of  the  flow  field.[19]  Principally,  LESLIE3D  autonomously  switches  on 
the  upwind  scheme  only  in  the  vicinity  of  flow  field  discontinuities  and  employs  the 
MacCormack  solver  in  smooth  regions  of  the  flow  field.  This  hybrid  methodology  performs 
quite  well  and  has  been  extensively  validated  for  the  Richtmyer- Meshkov  instability. [19]  It  has 
also  shown  a  great  deal  of  efficacy  for  simulating  shock-turbulence  interaction  problems  such  as 
the  reactive  blast  wave. 


3.3  Time  Integration  Scheme  and  Computer  Code  Structure 

LESLIE3D  utilizes  explicit  time  integration  to  capture  the  physics  associated  with  unsteady 
flow  fields.  Specifically,  we  find  that  explicit  time  integration  is  very  effective  for  accurately 
simulating  wave  propagation.  Although,  the  explicit  time  step  is  limited  by  the  Courant- 
Friedrich-Lewy  (CFL)  criteria,  it  presents  a  high  level  of  efficiency  and  computational  simplicity 
on  parallel  machines.  A  version  of  the  modified  Euler  method  is  applied  for  the  MacCormack, 
upwind  and  hybrid  schemes.  This  time  integration  method  (in  two  steps)  is  briefly  summarized 
for  one  space  dimension  as  follows. [20] 


jj*  n — 1  .  .  P./V;f 1 

ui  =ui  +  A“ ' (4 
Ax 


(+1/2 


?Num,n—'\  ' 
i-M2 


(49) 


Un=- 

1  2 


u: + u 


n- 1 


At 

Ax 


/  Nu  m,  * 

(-C+1/2 


(7  Nun/,*  \ 

C 1/2  ) 


(50) 


Equations  (49)  and  (50)  require  minor  adaptations  for  use  with  problems  cast  in  three 
dimensions.  Basically,  the  length  term  Av  roughly  corresponds  to  a  ratio  of  the  cell  volume  to  an 
average  surface  area  for  the  interface  normal  to  the  coordinate  indexed  by  i .  The  advective  terms 
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in  (2)  through  (5)  are  easily  discretized  for  structured,  hexahedral  finite  volumes  while  the 
viscous  terms  require  a  basic  generalized  coordinate  transformation.  [13]  The  mathematical 
procedures  are  relatively  easy  to  accomplish  and  are  quite  robust. 

3.4  Algebraic  Grid  Generation 

At  the  beginning  of  this  report,  we  mentioned  that  a  principal  goal  of  this  work  is  to  develop  a 
rapid  method  of  generating  grids  for  a  relatively  simple  body  geometry.  Our  prior  work  in  this 
area  required  a  great  deal  of  grid  generation  effort  (over  200  blocks).  The  workload  associated 
with  making  even  the  most  changes  to  this  grid  (such  as  adding  points)  is  daunting.  In  this  case, 
we  employ  a  single  block  grid  for  a  simpler  geometry  that  possesses  the  desired  lobed  design. 
Since  the  body  shown  in  Figure  1  still  possesses  some  symmetry,  the  elected  plan  of  attack 
entails  generating  a  series  of  two-dimensional  grids  parallel  to  the  longitudinal  body  center. 
These  two-dimensional  slices  can  easily  be  joined  laterally  to  create  the  three-dimensional  mesh. 
Any  of  the  longitudinal  slices  can  be  described  by  four  arcs:  one  along  the  body  surface,  the 
second  along  the  farfield  arc.  The  third  is  from  the  center  point  on  the  node  to  the  farfield,  and 
the  fourth  is  from  the  center  point  on  the  tail  to  the  farfield.  These  arcs  are  illustrated  in  Figure  2. 


x2(i) 


Figure  2.  Arrangement  of  arcs  for  a  two-dimensional  grid  slice 

By  establishing  a  mapping  between  from  arc  1  to  arc  2  and  then  from  arc  3  to  arc  4,  a  two- 
dimensional  mesh  may  be  generated.  Transfinite  Interpolation  (TFI)  is  a  commonly  applied 
algebraic  technique  used  to  perform  this  task.  [21]  TFI  is  a  very  fast  and  easily  programmable 
procedure.  Of  course,  it  has  its  disadvantages.  TFI,  in  its  original  form,  offers  little  control  over 
point  spacing  say,  adjacent  to  the  body  surface.  Also,  it  offers  no  direct  means  of  enforcing 
orthogonal  intersections  of  grid  lines.  One  may  recall  that  highly  non-orthogonal  meshes  can 
result  in  increased  numerical  errors  in  CFD  simulations.  These  errors  can  be  particularly 
prevalent  in  the  discretization  of  viscous  terms.  For  this  reason,  a  variant  of  TFI  involving 
Hermite  polynomials  is  utilized  for  this  work. 

In  its  first  incarnation,  TFI  provides  a  smooth  grid  (in  many  cases)  with  some  control  over 
point  spacing.  Unfortunately,  it  offered  no  control  over  the  slope  of  grid  lines  intersecting  the 
boundary.  This  problem  occurred  largely  because  the  blending  functions  employed  by  TFI 
allowed  only  the  specification  of  a  point  on  the  boundary. [21]  The  slope  at  this  point  could  not 
be  specified.  Later  on,  TFI  was  improved  by  incorporating  Hermite  polynomials  as  blending 
functions. [22]  These  polynomials  are  more  evolved  than  the  older  linear  blending  functions  and 
allow  a  value  for  the  slope  to  be  specified  at  the  endpoint  of  an  arc.  Although  offering  great 
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utility,  the  use  of  Hermite  polynomials  does  complicate  the  mapping  equations.  Referring  to 
Figure  2,  we  employ  the  four-boundary  method  developed  by  Shih  et  al.[ 22]  We  begin  by 
selecting  pairs  of  opposing  arcs  for  the  first  interpolation  step.  We  select  arcs  1  and  2  (“parallel” 
to  the  body  surface.  Let  points  along  arc  1  be  denoted  as  The  first 

interpolation  is  conducted  in  the  //  direction,  so  the  relevant  slopes  along  arc  1  are  designated  as 


dXAZ,0)  J  dYAZ, 0)  A11  r  J  ,  ,  A  J.  , 

and  — — — - .  All  of  these  values  must  be  supplied  by  the  user.  Accordingly, 


di] 


drj 


dX  (Z  1)  / 

analogous  information  must  be  provided  along  arc  2,  i.e.,  X2(^,1),  7’  and 

™%.  Let  two-dimensional  grid  points  in  the  field  for  step  1  be  designated  as  (jc1, y1); 
then  the  step  1  equations  are 


x1  (Z,  rj)  =  X  1  (£,  7)  /q  (77)  +  X2  (£,  rj)  h2  (7)  + 


dx(Z,rj  =  0) 
drj 


h 3(7)  + 


dx(Z,i)  =  X) 
drj 


Kin) 


.v1  (Z,  v)  =  (%>  if)  K  fi)  +  y2  (£>  if  h2  (u)  +  0)  h3  ff  +  ^  K  (u) 


(51) 


dr/ 


drj 


Functions  h  1  through  I14  are  the  Hermite  polynomials 

/q  =  27s  -3r/2  +1;  7i2  = -2r/3 +3r/2 
h3  =  r/3 -2r/2 +  r/ ;  7i4=r/3-r/2 


(52) 


Equations  (51)  contains  a  few  strange  quantities,  the  partial  derivatives  of  x  and  y.  These 
quantities  x  and  y  are  actually  point  coordinates  for  the  final  grid.  What  makes  these  equations 
strange  is  that  the  final  grid  coordinates  are  unknown  at  this  stage  of  the  procedure.  For  the  first 
two  or  three  readings  of  the  core  reference,  this  notation  proved  to  be  confusing.  What  one  must 
realize  is  that  the  partial  derivatives  in  (51)  are  strictly  defined  along  the  boundary  arcs  for  the 
final  grid.  Yet,  these  boundary  arcs  mathematically  coincide  (up  to  a  discrete  distribution  of  arc 
points)  with  arcs  (X^F^and  (X2,Y2).  Therefore,  the  notation  used  in  (51)  does  make  sense, 
and  the  properties  (curvature,  etc.)  of  the  curves 

X,:  (x(Z,i7  =  0),y(£, 7  =  0))  and  X2:  (x(Z,rj  =  1),y(^,7  =  1)) 

are  known  a  priori.  In  fact,  these  curves  provide  the  opportunity  to  control  the  slopes  of  grid 
lines  along  the  boundary.  We  may  specify  them  as  follows: 
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(53) 


o?i  dg  or/  dg 

_ft««  *§f«L. 

dr/  dg  dr/  dg 


Functions  (p^  2  are  used  to  control  orthogonality  of  the  grid  lines  and  also  to  help  prevent  grid 

adjacent  lines  from  crossing. [23]  Notice  that  the  derivatives  of  x  (with  respect  to  q)  depend  upon 
the  derivatives  of  y  (with  respect  to  g).  A  similar  statement  can  be  made  for  derivatives  of  y.  Note 
that  the  right  hand  sides  equation  (53)  are  calculable  since  point  coordinates  for  arcs  1  and  2  are 
known.  The  negative  signs  shown  in  (53)  are  essential  in  controlling  grid  line  orthogonality  at 
the  boundary.  More  particularly,  recall  that  the  normal  line  slope  to  a  point  on  an  arc  is  given  by 
the  negative  reciprocal  of  arc  slope  at  the  point. 

Equations  (52)  and  (53)  perform  an  accurate  and  well  controlled  interpolation  between  arcs  1 
and  2  in  the  q  direction.  However,  the  resulting  distribution  of  points  fails  to  match  arcs  3  and  4. 
A  second  interpolation  representing  the  E,  direction  is  warranted.  As  is  dictated  by  original  TFI 
theory,  the  second  step  interpolates  the  errors  generated  by  the  first  step. [21]  This  dictum  is 
enforced  here,  but  the  equations  differ  slightly  from  the  classic  TFI  equations.  The  goal  of  this 
step  is  to  determine  increments  (  Ax(^,  r/) ,  Ay(£,  77) )  such  that 

*(£,*7)  =  x1(^,/7)  +  Av(^,77); 

y(€,i 7)  =  y\Z,r/)  +  {±y(%,r/) 


where  (x(^,r/), y(^,r/) )  are  the  final  grid  points’  coordinates  (unless  one  elects  to  redistribute 
points  along  the  grid  lines).  These  increments  are  calculated  as  follows. 


A  x(£ ,  77) 

f 

+ 

V 


=  [A3(77)-x1(^  =  0,77)]/75(^) 

dx(%  =  0,r/)  dx\%  =  0,77)^ 


+ 

+ 


[A4(77)-x1(^  =  1,77)]/76(^) 

'dxt£  =  \,r/)  dx\%  =  \,rj)' 

l  d£  ^  J 


(55) 


Ay(£,77)  =  [Y3(ij)-y\4  =  0,r/)]h5(4)  +  [Y4(r/)-y\^  =  1,r/)]h6(^) 


dy(<^  =  0  ,r/) 

{  H 


^(^  =  0,77)^ 

^  J 


h7(Z) 


+ 


'dy{£  =  \,rj) 

l  ^ 


dy\%  =  1.7)^ 

^  J 


(56) 


where  the  Hermite  polynomials  are  expressed  as 


h5=2f  -3<f+1;  /i6=-2«f +3<f  ; 
h7=f-2f+^-  h8  =  f-f 


(57) 


The  partial  derivatives  of  x  and  y  set  along  the  boundary  in  (56)  and  (57)  are  defined  below. 
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(58) 


a^  =  a7)=  3(7) 


5£ 

54^  =  1, 7) 


dr/ 


84 


=  f,W®;  ^4^  =  -f.W 


dr] 
dX4(77) 


5£  "  5t;  *  di] 

Equations  (55)  and  (56)  also  contain  derivatives  of  the  stage  1  grid  coordinates  (jc1,  y1).  The 
formulas  for  these  derivatives  are  produced  by  differentiating  (51).  We  simply  list  the  results. 


fix1  (£  =  0,77) 

5£ 


fix1  (£  =  1,77) 

5£ 


=  ^i(7) 
+  773(77) 

=  *i(7) 
+  /i3(7) 


dy\%  =  0.77) 

5£ 


5y1(^  =  l?7) 
5£ 


=  ^1(7) 

+  /i3(77) 

=  *i(7) 
+  773(7) 


5x(<^  =  0, 77  = 

=  0) 

5£ 

52x(<^  =  0,77 

=  0) 

d^dr/ 

5x(£  =  1,77  = 

0) 

5£ 

52x(<^  =  1,77 

=  0) 

5<^577 

„  5y(4:  =  0,7  = 

=  0) 

5£ 

52y(£  =  0,77 

=  0) 

dd,dt] 

dy(£  =  1,7  = 

0) 

5£ 

52y(£  =  1,7 

=  0) 

+  /72(?7) 


+  /74(77) 


+  h2(rj) 


+  K  (7) 


+  /72(77) 


+  *4(7) 


5x(<^  =  0,7  =  1) 

5? 

52x(<^  =  0,77  =  1) 

d%dr] 

54^  =  1,77  =  1) 

5£ 

52x(^  =  1,?7  =  1) 
d^dr/ 

5y(^  =  0,77  =  1) 

5£ 

52.y(£  =  0,77  =  1) 


dt^dr/ 


+  /72(?7) 


+  *4(7) 


dE,di] 

5yQ?  =1,7  =  1) 

5£ 

52y(«^  =  1,77  =  1) 
d^dr/ 


(59) 


(60) 


(61) 


(62) 


With  the  use  of  the  equations  presented  above,  a  grid  is  easily  produced  for  a  two-dimensional 
slice.  Derivatives  along  the  boundary  arcs  occurring  in  (53)  and  (58)  can  be  either  analytically  or 
numerically  calculated.  Since  we  are  resolving  a  viscous  flow  field,  maintaining  control  over  the 
grid  point  spacing  normal  to  the  body  is  important.  For  this  reason,  points  along  constant  £,  grid 
lines  are  redistributed  after  the  two-dimensional  grid  has  been  generated.  We  have  elected  to 
redistribute  the  points  via  hyperbolic  sine  arc  length  transformation.  [21]  Given  the  prevalence  of 
this  transformation  in  modem  grid  generation,  the  details  are  omitted  from  this  report. 

The  final  step  in  this  process  is  to  create  the  three-dimensional  grid.  Two-dimensional  slices 
are  generated  at  each  azimuthal  station  around  the  circumference  of  the  body.  The  aerodynamic 
body  has  five  lobes,  so  the  shapes  of  any  two  adjacent  slices  differ  from  one  another.  Common 
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longitudinal  (i)  and  normal  ( j)  indices  are  used  to  generate  these  two-dimensional  mesh  slices, 
so  connecting  them  through  the  azimuthal  index  (k)  to  form  hexahedral  volumes  is  not  difficult. 
This  method  produces  an  adjustable  mesh  very  quickly  and  even  permits  grid  lines  spacing 
controls  for  the  mesh  interior.  Although  our  mesh  is  not  highly  orthogonal,  cells  adjacent  to  the 
body  do  have  nice  aspect  ratios. 

3.5  Computation  of  Turbulent  Statistics 

A  topic  of  great  interest  in  the  present  work  is  an  investigation  of  turbulent  flow  near  the 
aerodynamic  body  surface.  Localized  regions  of  elevated  subgrid  kinetic  energy  (a  primary 
indicator  variable  for  turbulence)  are  examined  with  statistical  analysis  methods.  We  consider 
three  such  analyses:  (i)  a  velocity  correlation  matrix,  (ii)  the  probability  distribution  function  for 
velocity  fluctuations  and  (iii)  the  spectrum  of  turbulent  kinetic  energy.  These  statistical 
investigations  are  not  exhaustive,  yet  they  do  render  insight  into  properties  of  the  turbulent  flow 
field  at  high  Mach  numbers.  To  perform  a  more  in-depth  study,  it  would  also  be  necessary  to 
analyze  the  wave  number  content  of  turbulent  kinetic  energy  for  two-point  correlations.  [24] 

The  velocity  correlation  matrix  shows,  on  the  average,  how  different  fluctuating  velocity 
components  correlate  with  one  another  in  time;  that  is  to  say,  are  the  fluctuations,  existing  at 
point  in  space,  “like”  one  another  in  the  way  they  change?  To  produce  the  matrix,  we  rely  on  the 
time  average  defined  as 


-|  'o+T 

u  —  —  J«(0  dt  (63) 

1 0 

Equation  (63)  is  interesting  because,  at  face  value,  it  seems  to  remove  all  time  dependency  from 
the  time  averaged  quantity.  In  fact,  this  interpretation  is  used  within  the  context  of  this  report. 
Still,  there  is  a  wider  interpretation  of  (63).  Since  the  time  average  relies  upon  parameters  to  and 
T,  to  is  the  point  in  time  where  the  average  is  anchored  while  T  can  be  thought  of  as  the  “width” 
of  the  averaging  window.  The  interpretation  is  more  clearly  illustrated  by  rewriting  (63)  with  the 
transformation  70=t0+T / 2 ,  i.e., 


1  io-r/2 

u(JQ)  =  —  ju(t)dt  (64) 

T  t0-T  12 

Velocity  fluctuations  are  computed  as  follows.  Let  «,  denote  the  ith  component  of  fluid 
velocity.  Then  the  average  of  this  component  is  denoted  as  if, ,  and  its  associated  fluctuation  is 
computed  as 


u'i=ui—ui  (65) 

Based  upon  (63),  the  velocity  (one -point)  correlation  matrix  is  defined  as 
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(66) 


u'  u';  =  u.  u ,  —  u  u . 

1  J  1  J  1  J 

where  the  properties  of  averaging  have  been  used.  [24]  Important  aspects  of  turbulent  physics  are 
contained  in  this  matrix.  For  instance,  the  turbulent  kinetic  energy  per  unit  mass  E  is  closely 
related  to  its  trace.  Observe  that 


E  =  ^Tr  [u'u'j]  =  +  u'2  +  u'22  )  (67) 

Equation  (67)  is  an  analogue  of  the  common  fluid  kinetic  energy.  The  off  diagonal  matrix  entries 
describe  how  changes  in  component  velocity  fluctuations  are  related.  For  example,  if  u\  u'2  >  0 . 
then  a  positive  change  in  u\  tends  to  coincide,  on  the  average,  with  a  positive  change  in  u2  .  On 

the  other  hand,  u\ u2  >0,  implies  the  changes  must  have  opposing  signs. [25]  These  correlations 
may  be  used  to  determine,  in  part,  the  level  homogeneity  or  isotropicity  possessed  by  turbulent 
fluctuations. 

Since  turbulent  fluctuations  are  inherently  random,  they  may  be  analyzed  in  accordance  with 
the  rules  of  probability.  Velocity  fluctuations  existing  at  a  point  in  space  are  therefore  distributed 
in  accordance  with  a  probability  distribution  function  (PDF).  Given  a  random  data  set,  we  may 
estimate  the  PDF  denoted  by  the  letter/.  The  simple  method  for  doing  so  is  by  histogram. [26] 
Suppose  that  the  data  consists  of  a  finite  number  N  of  bounded  velocity  fluctuations,  i.e.,  there 
exist  finite  numbers  u'  min  and  u\  m!lv  such  that 

(  ,  IlUil  l  ,  lllaA 


It. 


<  u'  <  u'. 


(68) 


for  all  N  samples.  Equation  (68)  delineates  a  domain  in  $R3  that  may  be  divided  into  a  number  of 
bins.  Each  bin  has  size  Aw'  =  ( u'j  ma  -  u'i  mn )  /  M t  where  the  M,-  are  the  numbers  of  bins 
established  by  the  analyst  along  each  coordinate.  In  this  report,  we  calculate  bivariate  joint  PDFs, 
e.g.,  /(«',«' ) .  The  mechanics  of  calculating  this  PDF  begin  by  setting  all  (M;)(M,j  bins  to  zero. 

Next  we  loop  over  all  N  fluctuations;  if  a  particular  pair  of  fluctuations  ( u' ,  w' )  falls  within  the  a 
bin  delineated  by  (68),  we  add  the  quantity 


A/  = 


1 

N  (Au')(Au'j) 


(69) 


to  the  bin  sum  for/.  Clearly,  for  samples  consisting  of  small  numbers  N  of  fluctuations,  many  bin 
sums  may  remain  zero,  so  /  may  be  discontinuous.  Still,  with  the  use  of  (68)  and  (69),  the 
cumulative  distribution  function  (CDF)  is  properly  normalized,  i.e., 


lim  F {u ' ,  u ' )  =  f  [  /(«',//  )  du'  du'  =  1 

>00  J  J  J  ‘  ' 


(70) 
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The  simple  PDF  is  adequate  for  studying  the  behavior  of  turbulent  velocity  correlations,  but  it  is 
not  very  smooth  (especially  for  small  data  sets).  A  smoother  PDF  may  be  constructed  by 
dispensing  with  the  idea  of  bins  and  summing  Gaussian  functions  for  each  pair  of  velocity 
fluctuations  in  the  data  set.  As  is  shown  below,  this  PDF  is  easily  written  in  the  form  of  an 
equation. 


/(M',v')  =  X^exp| 

n= 1  ’ 


I \u’-u'nf  (v*-p2V 
1  2of  2cr2  J 


(71) 


In  (71),  ( u'n  ,v')  is  a  pair  of  velocity  fluctuations  taken  from  the  data.  The  coordinate  widths  of 
the  Gaussian  are  taken  from  the  diagonal  entries  of  the  velocity  correlation  matrix  (66). 

Specifically,  <ju  =  ,  and  cjv  =  .  Equation  (71)  must  be  normalized  to  satisfy  the 

requirements  of  the  CDF  (70);  accordingly,  A  is  the  normalization  constant  and  can  be  shown  to 
be 


A  =  — ! -  (72) 

2/rau  av 

For  this  PDF,  selection  of  the  domain  of  integration  in  (70)  is  very  important.  During  numerical 
integration,  one  must  ensure  the  domain  is  large  enough  to  ensure  that  the  CDF  is 
computationally  normalized. 

The  final  tool  addressing  turbulent  statistics  employed  in  this  work  is  spectral  analysis. 
Turbulent  fluid  motions  are  characterized  by  a  wide  range  of  scales,  both  temporal  and  spatial. 
Turbulent  eddies  have  different  length  and  time  scales.  Spatial  scales  have  a  direct  spectral 
analogue  known  as  wave  number,  and  time  scales  have  a  similar  analogue  known  as  frequency. 
Both  of  these  analogues  lend  insight  to  the  structure  of  turbulence.  Here,  we  address  only  the 
time-frequency  spectrum  of  turbulent  kinetic  energy.  In  our  formulation,  turbulent  kinetic  energy 
E  is  defined  as  a  function  of  time  at  a  single  point  in  space,  i.e., 


E(t)  =  i(nf(r)  +  n'2(r)  +  n'2(r))  (73) 

The  frequency  spectrum  is  rendered  by  Fourier  transformation  of  (67). [24,25]  The  forward 
transformation  is  defined  as 


i  00 

«  |  |* 

E(a>)  =  —  E(t)exp(-ia>t)dt 
2  n  J 


while  the  inverse  transformation  is  written 


(74) 
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(75) 


oo 

E{t)=  J  E(o))exp(io)t)do) 

—oo 


where  co  is  angular  frequency  and  “i  ”  is  the  imaginary  unit.  Transforms  (74)  and  (75)  are 
“general”;  in  practice,  we  do  not  recognize  the  concept  of  negative  time  or  negative  frequency, 
so  the  lower  the  limits  of  the  transform  integrals  are  set  to  zero.  The  theoretical  transforms  are 
also  continuously  defined  unlike  the  discrete  data  set  addressed  later  in  this  report.  It  follows  that 
we  must  use  the  discrete  Fourier  transform  (DFT)  in  lieu  of  (74).  The  DFT  is  defined  for  N  data 
points  non-uniformly  distributed  in  time  as 


N- 1 


E(co)  =  —  YE(tn ) exp(— / cotn ) (fB+1  -tn),  co>  0 

C-71  n= 1 


(76) 


E(co )  is  a  complex  function  of  real  co,  so  to  express  the  frequency  content  of  E,  we  plot  the 
power  spectral  density  defined,  in  this  case  as 

Pe{(o)  =  4e{(o)*E{(o)  (77) 

a  real  valued  function  of  co.  DFTs  must  be  used  with  great  care,  especially  for  small  data  sets  due 
to  induced  numerical  errors  such  as  aliasing.  To  remedy,  in  part,  this  type  of  difficulty,  we 
observe  the  limitation  imposed  by  the  Nyquist  sampling  frequency. [27]  Noting  that  <x>=2nf , 
where/  is  regular  frequency  in  Hertz,  the  Nyquist  limitation  states  that  the  DFT  can  resolve  no 
frequencies  exceeding  fNyq  where 


■f Nyq 


2  At 


(78) 


where  Atirax  is  the  maximum  time  step  between  any  two  events  adjacent  in  the  time  series. 
Obviously,  from  (78)  there  is  inverse  relationship  between  the  time  step  between  data  points  and 
the  maximum  resolvable  frequency.  Secondly,  and  not  evidenced  by  (78)  are  relationships 
involving  the  minimal  resolvable  frequency  as  well  as  the  ability  to  resolve  two  adjacent 
frequencies.  This  type  of  resolution  say,  the  minimum  Af  is  inversely  proportional  to  the  length 
of  time  encompassed  by  the  entire  time  series.  Thus,  to  resolve  progressively  lower  frequencies, 
we  must  sample  the  time  series  for  more  time.  Moreover,  to  resolve  two  different  frequencies, 
we  must  sample  an  interval  of  time  that  covers  at  minimum  an  integer  number  of  complete  event 
cycles  for  both  frequencies.  Comparing  frequency  spectra  for  turbulent  kinetic  energy  at 
different  space  points  is  a  goal  of  this  report. 
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4  SETTING  UP  THE  SIMULATION 


This  section  of  the  report  describes  measures  taken  to  set  up  the  simulation.  We  endeavor  to 
include  enough  information  to  allow  others  to  conduct  their  own  numerical  studies  on  our 
geometry. 

4.1  Grid  Geometry 

A  picture  of  the  body  geometry  is  shown  in  Figure  1.  As  mentioned  in  Section  3.4,  the  grid  is 
constructed  from  two-dimensional  slices  connected  in  azimuth.  The  surface  of  the  body  is  the 
no-slip  surface  located  at  j  =  1.  The  overall  body  length  is  7  cm.  The  maximum  radius  of  the 
body  (rmax)  is  1.5  cm  while  the  minimum  body  radius  (rmin)  is  0.5  cm.  The  body  is  characterized 
by  five  lobes  arranged  in  azimuth  as  shown  in  Figure  3.  Because  of  the  lobes,  the  maximum 


Figure  3.  Azimuthal  plot  of  the  outer  contour  of  the  body  illustrating  the  lobe  positions. 

inner  radius  (contour  A1  in  Figure  2)  of  the  two-dimensional  slice  changes  at  each  azimuthal 
station.  The  variation  in  this  radius  is  dictated  by  the  following  equations. 

h  =h  +  Ah cos(5$) ,  0  <  6  <  2k 

h  =Hrnm+rnn)  (79) 

Ak  =Krmax  -rmin) 

The  front  or  nose  of  each  slice  is  comprised  of  a  disk  of  fixed  radius  rmm.  On  the  other  hand,  the 
rear  of  each  slice  is  comprised  of  a  disk  of  radius  rmax  only  when  6  =  2nn/5,  n  =  0,  1,  ...,  4.  In 
order  for  the  body  to  terminate  at  a  single  point,  slices  located  at  other  values  of  6  the  rear  slice 
contour  is  given  by  the  same  circular  arc  (radius  rmax)  multiplied  by  h.  This  multiplier  has  the 

effect  of  flattening  the  circular  arc.  Figure  4  illustrates  the  placement  of  three  of  the  X1  slice 
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contours;  note  that  the  vertices  for  the  front  and  rear  arcs  are  separated  by  five  centimeters.  It  is 
evident  that  arc  representing  the  rear  of  the  body  is  flattened  for  those  values  of  9  not  coinciding 

with  rmax.  The  final  portion  of  the  X1  contour  is  the  straight  line  segment  joining  the  forward  and 
rear  arcs.  Experience  has  shown  that  this  line  must  smoothly  intersect  (via  continuous 
derivative)  the  forward  and  rear  arcs  in  order  to  achieve  a  smooth  distribution  of  grid  lines. 


Figure  4.  Illustration  of  body  surface  contours  for  different  azimuth  angles  (slices).  Each  of  these  contours 
forms  the  Xt  contour  for  generating  the  grid  on  the  associated  azimuthal  slice. 

Smooth  intersections  may  be  accomplished  by  locating  two  points,  one  on  the  rear  side  of  the 
forward  arc  and  the  other  on  the  forward  side  of  the  rear  arc,  where  the  derivative  taken  with 
respect  to  the  longitudinal  coordinate  x  are  the  same.  By  including  more  points  along  this  contour 
the  intersection  becomes  progressively  smoother. 

The  far  field  grid  contour  delineated  by  points  X2  in  Figure  2  is  formed  by  a  simple  circular 
arc  with  a  radius  of  15  cm.  Contours  X3  and  X 4  are  straight  line  segments  that  are  common  to 

all  slices.  Note  also  that  the  overall  domain  has  a  spherical  shape.  These  four  contours  are  used 
as  inputs  by  the  Hermite  TFI  routine  described  in  Section  3.4.  The  grid  utilized  to  generate 
results  for  this  problem  has  301  points  in  the  longitudinal  (?)  direction  and  201  points  in  each  of 
the  normal  (/)  and  azimuthal  (k)  directions.  Minimum  spacing  at  the  body  surface  is  set  at  10  5 
meters  corresponding  to  a  y+  value  of  just  under  ten.  [24] 

4.2  Initial  Conditions 

As  mentioned  earlier,  the  aerodynamic  body  is  immersed  in  a  Mach  10  flow  field  at  sea  level. 
The  body  is  oriented  at  5°  pitch  and  5°  yaw.  The  simulation  is  started  by  setting  the  associated 
Cartesian  velocities  for  these  conditions  in  the  outer  spherical  shell  of  farfield  cells.  In  the  same 
cells,  subgrid  kinetic  energy  is  initialized  to  a  value  of  0.01  times  the  freestream  kinetic  energy. 
The  freestream  temperature  is  set  at  300°K.  Farfield  boundary  conditions  must  be  specified  for 
FESFIE3D.  This  process  is  performed  in  the  initial  conditions  by  checking  for  the  presence  of 

inflow  in  the  outer  cell  layer.  Suppose  the  V  is  the  local  velocity  and  n  is  the  outward  pointing 
normal  vector  for  a  farfield  cell.  The  boundary  condition  is  then  set  as 
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(80) 


Farfield  BC  = 


J  Inflow 
[Outflow 


V  »n<  0 

V  •  n >  0 


The  present  computational  study  uses  turbulent  mixing  and  is  non-reactive.  Since  it  is  non¬ 
reactive,  only  oxygen  and  nitrogen  gases  are  involved  in  the  problem.  This  problem  is  executed 
on  512  processors  (actually  cores).  The  average  time  step  is  about  1011  seconds,  and  1.383 
million  time  steps  are  executed.  Each  step  requires  about  0.7  seconds  of  wall  clock  time.  We 
allow  the  flow  to  settle  through  a  travel  distance  of  11.5  body  lengths  before  collecting  data. 


5  RESULTS 


The  lobed  body’s  Mach  6  flight  has  been  simulated  as  specified  in  the  Section  4.  Bulk 
properties  such  as  pressure,  temperature,  subgrid  kinetic  energy  and  vorticity  have  been  extracted 
from  the  numerical  solution  and  used  to  guide  the  extraction  of  statistical  information.  Time 
series  information  is  also  collected  regarding  oscillatory  forces  exerted  on  the  body  as  well  as 
turbulent  velocity  fluctuations  and  kinetic  energy. 
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Figure  5.  Contour  plots  of  vorticity  magnitude  set  on  a  vertically  oriented  plane  at  the  body  centerline.  Plot 


(a)  occurs  at  0.428  ms  solution  time  while  plot  (b)  occurs  at  0.439  ms. 


5.1  Surface  Vorticity 


Figure  5  shows  two  contour  plots  of  vorticity  magnitude  adjacent  to  the  body’s  nose  The  bow 
shock  is  clearly  visible  and  has  a  standoff  distance  of  0.678  mm  from  the  nose.  As  is  expected, 
high  levels  of  vorticity  are  generated  behind  the  curved  bow  shock  and  are  washed  downstream 
along  the  body  surface.  The  vorticity  field  undergoes  significant  change  over  the  0.011  ms  time 
interval  between  these  two  snapshots.  The  presence  of  strong  vorticity  fluctuations  in  this  region 
motivates  a  closer  examination  of  the  flow  region  containing  the  boundary  layer.  The  same 
vertical  plane  in  the  solution  is  examined  in  Figure  6  with  enhanced  magnification.  The  mesh  is 
overlain  on  these  plots  to  provide  a  means  of  geometric  reference.  These  plots  are  interesting 
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Figure  6.  Magnified  contour  plots  of  vorticity  magnitude  showing  eddy  organization  set  on 
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a  vertically 


oriented  plane  at  the  body  centerline.  Plot  (a)  occurs  at  0.428  ms  solution  time  while  plot  (b)  occurs  at  0.439 


ms. 


Longitude  Station  (downstream  from  nose) 

Figure  7.  Vorticity  magnitude  plotted  for  the  body  adjacent  cell  at  seven  different  solution  times.  The 
horizontal  axis  relates  the  number  of  grid  cells  downstream  from  the  node  along  the  body  surface. 


since  we  can  see  vorticity  generated  at  the  shock  surface  (as  expected),  yet  it  is  even  more 
compelling  to  observe  the  flow  in  the  region  nearest  to  the  body  surface.  Without  losing  too 
much  generality,  we  can  regard  the  regions  inside  of  closed  contours  in  6(a)  and  6(b)  as  eddies. 
Nearest  the  wall,  these  eddies  seem  to  propagate  upstream  (to  figure  left)  and  then  move 
downstream  along  the  body.  This  motion  is  certainly  admissible  since  the  flow  behind  the  curved 
shock  is  subsonic.  Due  to  the  coupling  of  vorticity  and  pressure  fluctuations,  we  expect  that  a 
form  of  feedback  exists  in  this  region  of  the  flow,  so  it  becomes  worthwhile  to  assess  the  strength 
of  vorticity  fluctuations  in  the  boundary  layer.  These  fluctuations  may  drive  the  production  of 
turbulence.  Figure  7  provides  an  answer  to  this  question  albeit  over  a  fairly  limited  time  interval. 
This  figure  plots  vorticity  magnitude  in  the  body  adjacent  cell  over  a  span  of  0.141  ms.  Each 
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curve  plots  vorticity  magnitude  at  each  station  along  the  body  surface  from  the  nose.  It  is  evident 
that  boundary  layer  vorticity  is  very  unsteady  in  time  and  drives  the  production  of  subgrid  kinetic 
energy  through  contributions  to  total  strain  rate  in  equation  (24).  As  a  result,  turbulence  is 
produced  in  the  boundary  layer  and  continues  to  affect  the  downstream  flow  field. 

5.2  Surface  Pressure  and  Temperature 
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Figure  8.  Plots  of  thermodynamic  pressure  at  the  body  surface  for  (a)  0.428  ms  and  (b)  0.439  ms. 

Hypersonic  flow  fields  are  characterized  by  escalating  pressure  and  temperature  at  the  surface 
of  an  immersed  body.  Pressure  rises  sharply  at  stagnation  points  in  the  flow  field  while  the 
temperature  rise  is  more  widespread  as  the  flow  decelerates  under  the  force  of  skin  friction. 
Figure  8  contains  plots  of  pressure  at  the  body  surface  at  0.428  ms  and  0.439  ms.  Given  the 
relatively  simple  geometric  configuration  for  the  body,  the  pressure  field  is  unremarkable  except 
for  pressure  fluctuations  on  the  nose.  Due  the  action  of  turbulent  eddies  at  the  surface,  pressure 
contours  distort  and  change  periodically  in  time.  This  distortion  is  also  noted  on  the  lower  left 
“windward”  lobe  of  the  configuration  where  the  oncoming  hypersonic  flow  strikes  head  on.  The 
numerical  solutions  show  that  under  time  averaging,  pressure  fluctuations  are  strongly  correlated 
with  fluctuations  in  the  velocity  components.  Table  1  contains  estimates  for  the  pressure-velocity 
correlations  calculated  at  three  points  near  the  tip  of  the  nose. 

Table  1.  Pressure-Velocity  correlation  estimates  for  sampling  points  at  the  body  surface  near  the  nose  tip. 


Top 

Leeward 

Windward 

<  Pu  > 

1.71  x  10'4 

9.65  x  10'4 

2.09  x  10'3 

<Pv  > 

3.55  x  1CV 

1.36  x  10‘" 

-1.27  x  10"2 

<Pw> 

3.99  x  lO'3 

1.57  x  KT3 

-7.93  x  10'J 

Vertical  velocity  v  shows  the  strongest  correlation  with  pressure  fluctuations  at  sampling  point 
on  the  “windward”  side  of  the  nose  relatively  close  to  the  centerline. 

For  this  body  configuration,  both  surface  temperature  and  temperature  gradient  have  been 
examined.  Temperature  is  of  particular  interest  due  to  its  role  in  inducing  chemical  reactions 

25 

Distribution  A 


between  oxygen  and  nitrogen  in  the  surrounding  air.  Although  we  have  not  yet  attempted  to 
simulate  the  attendant  chemical  reactions,  we  can  determine  as  to  whether  or  not  temperature  is 
high  enough  to  support  ionization  and  dissociation  reactions.  Body  surface  temperature  based 
upon  the  adiabatic  wall  assumption  is  shown  in  Figure  9.  These  plots  suggest  that  a  significant 


Figure  9.  Plots  of  surface  temperature  (°K)  based  upon  the  adiabatic  wall  assumption  set  at  solution  times  (a) 
0.428  ms  and  (b)  0.439  ms. 


temperature  gradient  exists  in  the  lobed  region  of  the  body.  We  can  also  see  evidence  of 
temperature  fluctuations  over  the  body  surface  (analogous  to  the  pressure  fluctuations  noted 
earlier).  Figure  10  contains  plots  of  the  temperature  gradient’s  magnitude  calculated  at  the  body 
surface.  These  plots  do  show  extremely  high  temperature  gradients  at  the  body  surface  particu- 
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Figure  10.  Plots  of  surface  temperature  gradient  magnitude  (°K/m)  based  upon  the  adiabatic  wall  assumption 
set  at  solution  times  (a)  0.428  ms  and  (b)  0.439  ms. 


larly  on  the  windward  lobes  and  even  in  the  shielded  region  between  the  leeward  lobes.  In  the 
hottest  regions,  temperatures  exceed  2000°K.  As  a  result,  chemical  reactions  are  permissible  over 
a  large  fraction  of  the  frontal  body  area. 
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5.3  Turbulence  Statistics 


In  keeping  with  previous  remarks,  a  primary  focus  of  this  numerical  study  is  to  examine  the 
evolution  of  turbulent  flow  at  hypersonic  Mach  number  for  an  immersed  body  of  relatively 
simple  construction.  For  LDKM,  subgrid  kinetic  energy  drives  the  production  of  subgrid  stresses 
and  therefore  turbulence.  We  begin  by  examining  plots  of  subgrid  kinetic  energy  shown  in 
Figure  11.  Two  points  on  the  body  surface  are  selected  in  Region  1  for  statistical  sampling,  and 


Figure  11.  Plots  of  subgrid  kinetic  energy  on  the  body  surface  at  0.428  ms  (a)  frontal  view  and  (b)  oblique 
view.  Points  are  selected  for  statistical  sampling  from  region  1  and  from  region  2  on  the  lobe  surfaces. 

likewise,  two  points  are  selected  in  Region  2  where  subgrid  kinetic  energy  is  substantially 
elevated.  At  each  sampling  point,  time  series  are  computed  for  velocity  component  fluctuations 
as  well  as  for  turbulent  kinetic  energy.  Velocity  correlation  matrices  are  also  computed  for  the 
velocity  fluctuations.  In  concert  with  the  time  series,  this  information  is  used  to  compute 
probability  distribution  functions  (PDFs)  for  velocity  fluctuations  at  each  sampling  point.  These 
PDFs  are  presented  and  discussed  below. 

The  first  sampling  point  is  chosen  at  (i,k)  =  (200,76)  in  Region  1.  Both  simple  and  Gaussian 
PDFs  for  the  nv-correlation  have  been  obtained  at  this  point  and  are  presented  in  Figure  12. 
PDFs  for  the  uw  and  vw  correlations  are  presented  in  Figures  13  and  14,  respectively.  The  simple 
PDF,  indicated  by  (a)  in  each  figure  very  much  resembles  a  scatter  plot  formed  for  the  velocity 
fluctuations  and,  in  itself,  allows  the  analyst  to  determine  the  sense  of  the  velocity  correlation. 
On  the  other  hand,  the  Gaussian  PDF  serves  as  better  indicator  of  symmetry  or  of  skewness 
associated  with  the  distribution  function.  These  three  figures  are  interesting  because  all  of  the 
correlations  are  positive.  An  increase  (or  decrease)  in  one  velocity  component  tends  to  indicate  a 
respective  increase  (or  decrease)  in  the  other  two  velocity  components.  Also,  the  PDFs  indicate  a 
small  skewness  in  the  distributions.  The  mean  is  shifted  to  the  “right”  in  each  case.  Perhaps  due 
to  limitations  in  the  amount  of  data,  the  distributions  are  not  fully  symmetric.  An  increase  in  the 
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Figure  12.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uv  correlation  at  Region  1  sampling  point  (200,76). 


Figure  13.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uw  correlation  at  Region  1  sampling  point  (200,76). 
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Figure  14.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  vw  correlation  at  Region  1  sampling  point  (200,76) 


length  of  the  time  series  may  change  this  result. 


PDFs  for  the  Region  1  sampling  point  (i,k)  =  (210,75)  are  provided  in  Figures  15,  16  and  17. 
These  figures  contain  PDFs  for  the  uv,  uw  and  vw  velocity  correlations,  respectively.  These  plots 
are  interesting  because,  even  though  the  sampling  point  is  geometrically  close  to  (200,76),  the 
uw  and  vw  correlations  are  negative.  This  result  is  the  opposite  of  that  encountered  at  the 
neighboring  sample  point.  Secondly,  these  two  correlations  exhibit  a  high  level  of  skewness; 
their  PDFs  are  noticeably  stretched.  These  results  indicate  that  the  eddies  with  drastically 
opposite  velocity  arrangements  occur  very  close  together  on  the  body  surface. 
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Figure  15.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uv  correlation  at  Region  1  sampling  point  (210,75). 


Figure  16.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uw  correlation  at  Region  1  sampling  point  (210,75). 
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Figure  17.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  vw  correlation  at  Region  1  sampling  point  (210,75). 


In  Figures  18,  19  and  20,  the  velocity  correlation  PDFs  (retaining  the  same  order)  are  shown 
for  sampling  point  (i,k)  equal  to  (215,166).  This  sampling  point  resides  in  Region  2  as  indicated 
in  Figure  11.  The  uv  and  uw  PDFs  shown  in  Figures  18  and  19  indicate  correlated  velocity 
components,  but  corresponding  correlation  matrix  entries  are  somewhat  low.  Still,  the  u  and  v 
fluctuations  are  positively  correlated  at  this  point.  The  corresponding  simple  PDFs  show  a 
significant  amount  of  scatter  in  the  solution  as  evidence  of  poorer  correlation.  On  the  other  hand, 
the  vw  correlation  is  very  positive  and  clearly  reflects  the  order  exhibited  by  the  simple  PDF. 
Each  of  these  PDFs  exhibits  a  skewed  distribution  tending  to  favor  a  decrease- v  —  decrease- w 
tail  in  the  PDF. 
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Figure  18.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uv  correlation  at  Region  2  sampling  point  (215,166). 


Figure  19.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uw  correlation  at  Region  2  sampling  point  (215,166). 
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In  the  way  of  clarification,  the  PDFs  shown  in  the  figures  in  this  section  of  the  report  are 
graphed  in  terms  of  equi-probability  contours.  This  statement  applies  to  both  the  simple  and 
Gaussian  PDFs.  Figures  21,23  and  23  exhibit  the  PDFs  for  sampling  point  (i,k)  set  at  (230,168), 
a  point  in  relatively  close  proximity  to  the  preceding  point  (215,166)  also  in  Region  2.  One  can 
see  that  the  v  velocity  fluctuations  correlate  negatively  with  the  corresponding  u  and  w 
components  while  the  u  and  w  positively  correlate.  It  is  also  evident  that  the  correlation  matrix 
components  are  quite  large  representing  large  covariance  values.  The  resulting  PDFs  are  spread 
widely  over  the  velocity  planes  and  have  reduced  point  magnitudes.  Each  of  these  PDFs  is 
skewed  with  a  significant  tail  leading  away  from  the  mean.  Yet,  there  is  less  scatter  among  the 
fluctuating  velocity  components.  Some  of  the  skewness  may  be  due  to  the  small  sample  size 
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Figure  21.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uv  correlation  at  Region  2  sampling  point  (230,168). 
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Figure  22.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  uw  correlation  at  Region  2  sampling  point  (230,168). 


Figure  23.  Simple  (a)  and  Gaussian  (b)  PDFs  for  the  vw  correlation  at  Region  2  sampling  point  (230,168). 

used  for  this  analysis.  Overall,  it  is  interesting  to  see  how  different  the  covariance  matrix  entries 
are  for  two  closely  spaced  points  in  the  same  region.  Examples  are  the  uw  correlations  for 
Region  1  and  both  the  uv  and  vw  correlations  calculated  in  Region  2.  At  each  sampling  point,  the 
velocity  fluctuations  exhibit  randomness  albeit  with  components  with  widely  differing 
magnitudes  in  certain  cases.  The  covariance  matrix  entries  illustrate  this  concept  quite  well  and 
prove  very  useful  in  computing  the  Gaussian-based  PDF. 
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5.4  Turbulence  Spectra 

Turbulence  involves  a  wide  range  of  scales  in  fluid  motion  as  energy  cascades  from  the  scales 
of  the  mean  (exterior)  flow  field  down  to  the  smallest  scales  where  fluid  kinetic  energy  is 
dissipated  in  the  form  of  heat.  [24]  Thus  far,  we  have  conveyed  no  direct  information  that 
indicates  the  scale(s)  associated  with  velocity  fluctuations.  In  this  section,  we  provide  evidence 
of  the  scales  involved  through  the  inclusion  of  time-based  spectra.  Here  we  examine  spectra  for 
velocity  component  fluctuations  and  of  turbulent  kinetic  energy  (67).  Power  spectral  densities 
are  calculated  via  discrete  Fourier  transformation  (Section  3.5)  to  convey  the  frequency  content 
of  fluctuating  quantities.  The  spectra  are  calculated  at  the  same  sampling  points  used  in  the 
preceding  section.  Consider  sampling  point  (i,k)  =  (200,76);  its  power  spectral  density  is 
provided  in  Figure  24.  Although  the  spectrum  is  inactive  for  frequencies  below  200,000  rad/s, 


Angular  Frequency  (rad/s) 

Figure  24.  Power  spectral  densities  for  all  three  fluctuating  velocity  components  and  for  turbulent  kinetic 
energy  calculated  at  sampling  point  (200,76)  in  Region  1. 

but  there  is  a  great  deal  of  activity  in  the  high  frequency  section.  In  its  entirety,  the  spectrum 
shown  occurs  below  the  Nyquist  frequency.  The  lack  of  activity  in  the  low  frequency  section  is 
likely  due  to  the  short  sampling  interval  of  0.049  ms.  A  significant  increase  in  the  sampling 
interval  is  needed  not  only  to  capture  lower  frequency  eddies  but  also  to  sharpen  or  better  resolve 
the  high  frequency  peaks  in  the  spectrum.  Figure  24  is  clear  evidence  of  a  great  deal  of  small 
scale  motion.  This  behavior  is  also  reflected  at  sampling  point  (210,75).  Its  spectrum,  again 
computed  in  Region  1,  is  shown  in  Figure  25.  The  notation  F(«),  F(v),  F(w)  and  F(TKE) 
represents  the  power  spectral  densities  for  the  three  fluctuating  velocity  components  and  for 
turbulent  kinetic  energy,  respectively.  Although  this  spectrum  also  has  low  activity  at  lower 
frequencies,  the  high  frequency  bands  are  even  more  energetic.  This  behavior  is  quite 
noteworthy  because  the  two  sampling  points  in  question  are  in  close  proximity  to  one  another. 
For  further  information,  we  consider  sampling  points  in  Region  2.  The  spectra  for  sampling  point 
(215,166)  are  presented  in  Figure  26.  Again,  the  high  frequency  segment  of  the  spectrum  is  very 
active,  but  the  low  frequency  segment  is  not  developed  well  due  to  the  short  sampling  interval. 
The  smooth  curves  in  this  section  are  not  characteristic  of  fully  developed  spectra.  A  jagged 
spectrum  is  expected  in  the  low  frequency  range,  but  a  greatly  expanded  sampling  interval  is 
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needed  in  order  to  obtain  it.  Similar  spectral  behavior  is  exhibited  in  Figure  27,  computed  for 
sampling  point  (230,168)  in  Region  2. 


Figure  25.  Power  spectral  densities  for  all  three  fluctuating  velocity  components  and  for  turbulent  kinetic 
energy  calculated  at  sampling  point  (210,75)  in  Region  1. 


Angular  Frequency  (rad/s) 

Figure  26.  Power  spectral  densities  for  all  three  fluctuating  velocity  components  and  for  turbulent  kinetic 
energy  calculated  at  sampling  point  (215,166)  in  Region  2. 

5.5  Aerodynamic  Forces 

From  LESLIE3D’s  numerical  solutions,  we  have  computed  the  aerodynamic  force 
components  Fx,  Fy  and  Fz  exerted  on  the  immersed  body.  This  information  is  also  amenable  to 
time  series  analysis;  in  particular,  frequencies  of  oscillation  can  be  extracted  from  the  force  time 
history.  The  turbulent,  viscous  flow  field  has  a  pronounced  effect  on  the  aerodynamics  of  the 
body  since  tractions  (skin  friction)  and  separated  flow  (at  higher  pitch/yaw  angles)  can  strongly 
alter  the  forces  exerted  on  the  body.  The  time  history  for  Fx  is  shown  in  Figure  28.  This  plot 
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Figure  27.  Power  spectral  densities  for  all  three  fluctuating  velocity  components  and  for  turbulent  kinetic 
energy  calculated  at  sampling  point  (230,168)  in  Region  2. 


Time  (ms) 

Figure  28.  Time  history  plot  of  the  jc-component  Fx  of  aerodynamic  force  exerted  on  the  body. 

a  great  deal  of  high  frequency  activity.  Yet,  there  is  evidence  of  some  low  frequency  activity  that 
may  not  be  fully  resolved  in  the  relatively  small  sampling  window  applied  here.  Some  flow 
separation,  not  yet  revealed,  may  cause  larger  excursions  in  Fx  due  to  vortex  shedding.  Plots  of 
the  time  histories  Fy  and  F-  are  shown  in  Figure  29.  In  the  y  and  z  directions,  the  force  excursions 
are  not  as  large  as  in  the  x  direction,  but  the  high  frequency  content  is  still  present.  A  low 
frequency  oscillation  may  also  be  present.  As  before,  our  sampling  interval  is  too  short  to  fully 
reveal  this  dynamic.  Power  spectral  densities  have  been  extracted  from  the  force  component  time 
histories.  Plots  of  frequency  content  for  these  force  components  are  shown  in  Figure  30.  In  spite 
of  the  small  sampling  window,  the  discrete  Fourier  transformation  does  resolve  some  frequency 
content.  A  fundamental  frequency  of  oscillation  is  revealed  at  82,000  rad/s.  A  number  of  higher 
frequency  peaks  are  also  captured  and  not  necessarily  at  integer  multiples  of  the  fundamental. 
For  instance,  a  secondary  peak  is  located  at  235,000  rad/s  and  a  third  at  389,000  rad/s.  The  high 
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Figure  29.  Time  history  plots  of  the  x  and  v-components  Fx  and  Fy  of  aerodynamic  force  exerted  on  the  body. 


Figure  30.  Power  spectral  density  plots  for  the  three  aerodynamic  force  components. 


frequency  section  of  the  spectrum  is  very  active  for  all  three  components  indicating  the  presence 
of  many  oscillatory  frequency  all  below  the  Nyquist  sampling  frequency.  This  information 
stands  as  a  reasonable  indicator  for  the  presence  of  a  turbulent  flow  around  the  body. 

6  CONCLUSIONS 

This  report  has  addressed  some  aspects  of  the  hypersonic  aerodynamics  of  lobed  body.  In 
particular,  we  are  concerned  with  evolution  of  turbulent  flow  around  the  body  with  a  distribution 
of  subgrid  kinetic  energy  defined  at  the  far  field  boundary  of  the  computational  domain.  The 
simulation  is  conducted  at  Mach  6,  sea  level  flight  conditions,  with  pitch  and  yaw  angles  fixed  at 
five  degrees.  LESLIE3D,  a  multiphase  physics  computer  code  developed  by  Suresh  Menon  at 
Georgia  Tech,  is  used  to  perform  the  simulation  on  512  processors.  The  algorithms  employed  by 

35 

Distribution  A 


LESLIE3D  include  the  locally  dynamic  subgrid  kinetic  energy  model  with  fixed  coefficients.  A 
number  of  results  are  shown  in  this  report  such  as  body  surface  pressure,  temperature  and 
subgrid  kinetic  energy.  The  standoff  distance  for  the  curved  bow  shock  has  also  been  included, 
and  a  study  of  vorticity  evolution  behind  the  shock  and  in  the  boundary  layer  is  discussed. 
Statistics  for  the  turbulent  flow  field  have  been  calculated,  namely,  the  velocity  covariance 
matrix.  Probability  distribution  functions  have  been  computed  for  velocity  component 
fluctuations  and  turbulent  kinetic  energy.  Time-based  spectra  are  also  provided  for  the  same 
quantities.  Aerodynamic  force  time  histories  are  also  provided,  and  the  attendant  frequency 
content  has  been  computed. 

The  flow  field  is  characterized  by  turbulent  eddies,  particularly  in  regions  on  the  lobes  where 
subgrid  kinetic  energy  is  very  high.  The  time  sampling  interval  of  0.049  ms  is  rather  short,  so  the 
spectra  are  not  fully  developed,  particularly  at  lower  frequencies.  That  is  to  say,  to  resolve  a 
frequency  of  100  rad/s,  a  sampling  time  on  the  order  of  1.6  ms  is  required.  Clearly,  more 
simulation  time  is  required  to  resolve  frequencies  this  low.  Still,  the  velocity  field  contains  a 
great  deal  of  random  character,  and  the  high  frequency  spectral  band  is  very  active.  Fluctuating 
quantities  like  velocities,  pressure  and  aerodynamic  force  demonstrate  a  high  degree  of 
randomness  as  expected  for  a  turbulent  flow  field.  This  behavior  is  reflected  in  the  results  of 
spectral  analyses  performed  on  the  velocity  field  and  aerodynamic  force  components.  Each 
spectrum,  although  limited  by  a  relatively  short  sampling  window,  is  characterized  by  the 
presence  of  many  high  frequency  components  all  captured  below  the  Nyquist  frequency.  Further 
spectral  study  of  this  flow  field  is  warranted  for  a  substantially  wider  sampling  window. 
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